function [p_values_tests, significance_tests, rel_impulse_vector, F_test_all, incon_impulse_vector] = zFC_IdentifExtInstruments2(Resid,Instruments)

% given residuals from a reduced form var and some instruments, runs
% validity tests and computes the column of the B matrix, normalized to 1
% on the first entry. 

    % Make sure variables enter as column vectors
    if size(Resid,1) < size(Resid,2) 
        Resid = Resid'; 
    end
    if size(Instruments,1) < size(Instruments,2) 
        Instruments = Instruments'; 
    end

    % Make sure they have the same length
    assert(size(Resid,1) == size(Instruments,1))

    % Create operational variables
    obs           = size(Resid,1);
    n_variab      = size(Resid,2);
    n_instruments = size(Instruments,2);

    % Tests on validity of instruments
    p_values_tests     = NaN*ones(n_variab,4);
    significance_tests = NaN*ones(n_variab,4);
    F_test_all         = NaN*ones(n_variab,1);

    for i = 1:n_variab;
        yyy      = Resid(:,i);
        xxx_rest = ones(obs,1);
        xxx = [ones(obs,1), Instruments];

        [F_test, Autocorr_test, White_test]  = zFC_Ftest_Autocorr_Homosk(yyy,xxx,xxx_rest);
        F_test_all(i)       = F_test(1);
        p_values_tests(i,1) = F_test(2);
        p_values_tests(i,2) = Autocorr_test(1,2);
        p_values_tests(i,3) = Autocorr_test(2,2);
        p_values_tests(i,4) = White_test(2);


        significance_tests(i,1) = zFC_significance_stars(F_test(2));
        significance_tests(i,2) = zFC_significance_stars(Autocorr_test(1,2));
        significance_tests(i,3) = zFC_significance_stars(Autocorr_test(2,2));
        significance_tests(i,4) = zFC_significance_stars(White_test(2));
    end

    % Estimate relative impulse vector
    rel_impulse_vector    = NaN*ones(n_variab,1);
    rel_impulse_vector(1) = 1;

    yyy      = Resid(:,1);
    xxx      = [ones(obs,1), Instruments];
    beta_hat = (xxx'*xxx)^(-1)*xxx'*yyy;
    fitted   = xxx*beta_hat;

    for j = 2:n_variab
        yyy = Resid(:,j);
        xxx = [ones(obs,1), fitted];
        gamma_hat = (xxx'*xxx)^(-1)*xxx'*yyy;

        rel_impulse_vector(j) = gamma_hat(2);
    end


    yyy      = Resid;
    xxx      = [ones(obs,1), Instruments];
    beta_hat = (xxx'*xxx)^(-1)*xxx'*yyy;
    incon_impulse_vector = beta_hat(2,:)';

end